import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import contextily as ctx
import json
from matplotlib import rcParams
# 设置中文字体（SimHei 黑体）
rcParams['font.sans-serif'] = ['SimHei']  # 黑体
rcParams['axes.unicode_minus'] = False  # 解决负号问题

# 读取第1种格式的站点数据（从txt文件）
with open('SURF_CHN_MUL_MIN_202301010000.txt', 'r', encoding='utf-8') as file:
    # 跳过第一行标题
    next(file)
    # 读取数据
    data1 = [line.strip().split('\t') for line in file]

# 转换为DataFrame
df1 = pd.DataFrame(data1, columns=[
    "Station_Name", "Station_Id_C", "Station_Id_d", "Station_levl", "Station_Type", "Lat", "Lon", "Alti", 
    "Country", "Province", "City", "Cnty", "Town", "Admin_Code_CHN", "Year", "Mon", "Day", "Hour", "Min", 
    "PRS_Sensor_Alti", "PRS", "TEM", "RHU", "WIN_D_Avg_1mi", "WIN_S_Avg_1mi", "GST", "GST_5cm", "GST_10cm", 
    "GST_15cm", "GST_20cm", "GST_40cm", "LGST"
])

# 将Lat和Lon转换为数值类型
df1['Lat'] = df1['Lat'].astype(float)
df1['Lon'] = df1['Lon'].astype(float)

# 转换为GeoDataFrame
gdf1 = gpd.GeoDataFrame(df1, geometry=gpd.points_from_xy(df1.Lon, df1.Lat), crs="EPSG:4326")

# 读取第2种格式的站点数据（从json文件）
with open('tw.json', 'r', encoding='utf-8') as file:
    data2 = json.load(file)['RECORDS']

# 转换为DataFrame
df2 = pd.DataFrame(data2)

# 将lon和lat转换为数值类型
df2['lon'] = df2['lon'].astype(float)
df2['lat'] = df2['lat'].astype(float)

# 转换为GeoDataFrame
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2.lon, df2.lat), crs="EPSG:4326")

# 读取中国省边界数据
china_provinces = gpd.read_file('china.geojson')

# 绘制地图
fig, ax = plt.subplots(figsize=(12, 8))

# 绘制中国省边界
china_provinces.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=0.5)

# 绘制第1种站点数据
gdf1.plot(ax=ax, color='red', marker='o', markersize=1, label='站点1')

# 绘制第2种站点数据
gdf2.plot(ax=ax, color='red', marker='o', markersize=1, label='站点2')

# 添加底图
# ctx.add_basemap(ax, crs=gdf1.crs.to_string(), source=ctx.providers.Stamen.TonerLite)

# # 设置坐标轴范围
ax.set_xlim([108, 123])  # 根据经度范围设置
ax.set_ylim([18, 29])    # 根据纬度范围设置
ax.set_xlabel('经度', fontsize=12)
ax.set_ylabel('纬度', fontsize=12)

# 添加图例
# ax.legend()

# 设置标题
plt.title('气象站点分布图')

# 显示地图
plt.grid(True)
plt.show()